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ABSTRACT 

We use a set of cosmological N-body simulations to investigate the structural shape 
of galaxy-sized cold dark matter (CDM) halos. Unlike most previous work on the 
subject — which dealt with shapes as measured by the inertia tensor — we focus here 
on the shape of the gravitational potential, a quantity more directly relevant to com- 
parison with observational probes. A further advantage is that the potential is less 
sensitive to the effects of substructure and, as a consequence, the isopotential surfaces 
are typically smooth and well approximated by concentric ellipsoids. Our main result 
is that the asphericity of the potential increases rapidly towards the centre of the halo. 
The radial trend is more pronounced than expected from constant flattening in the 
mass distribution, and reflects a strong tendency for dark matter halos to become 
increasingly asphcrical inwards. Near the centre the halo potential is approximately 
prolate on average ((c/a)o = 0.72 ± 0.04, (b/a) — 0.78 ± 0.08), but it becomes less 
axisymmetric and more spherical in the outer regions. The principal axes of the isopo- 
tential surfaces remain well aligned, and in most halos the angular momentum tends 
to be parallel to the minor axis and perpendicular to the major axis. This suggests 
that galactic disks may form in a plane where the potential is elliptical and where 
its ellipticity varies rapidly with radius. This can result in significant deviations from 
circular motion in systems such as low surface brightness galaxies (LSBs) , even for rel- 
atively minor deviations from circular symmetry. Simulated long-slit rotation curves 
often appear similar to those of LSBs cited as evidence for constant density "cores" . 
This suggests that taking into account the 3D shape of the dark mass distribution 
might help to reconcile such evidence with the cuspy mass profile of CDM halos. 

Key words: cosmology: theory - cosmology: dark matter - galaxies: formation - 
galaxies: spiral - galaxies: kinematics and dynamics 



1 INTRODUCTION 



The fl at rotation curves of d isk g alaxies (|Rubin fe Fordl 
Il97d : iRoberts fc Whitehurstl 1 19751) . together with the 
anomalously h igh velocity dispersion of galaxies in clusters 
(Zwicky 1933), constituted the first compelling evidence for 
the existence of a massive dark matter component in the 
Universe. Since the interpretation of these data depends on 
the structure of dark matter halos, many theoretical studies 
have focused on the radial distribution of m ass within viri- 
alized dark matter halos, both analytically ([Cunn fc Gottl 
ll972l : lFillmore fc Goldreichlll984l:lHoffma,n fc Shaham|l985h 
and numerically llFrenk et al.lll985l.ll988 ; Quinn et all 19861 : 
iDubinski fc Carlberd Il99ll : I Crone et all Il994r t. This effort 
culminated in the realization that, in the prevailing cold 
dark matter (CDM) paradigm, the spherically averaged 



mass profile of dark halos is accurately described by scaling 
an ap proximately "universal" profile (|Navarro et al.l Il996l 
1 1991 NFW). The fitting formula proposed by NFW has 
enabled a simple means of testing the theoretical expec- 
tations against observations. Overall, these studies have 
shown that the mass profile of cold dark matter halos is 
broadly consistent with observations of X-ray clusters, grav- 
itation al lensing, satellite dynamic s, and galaxy group prop- 
erties dPointecouteau et al.l [20051; IComerford et al. I 120061 ; 



iPrada et al.ll2003l : iMandelbaum et ; aUl2006h . although prob- 
lems remain, especially on the scale of galaxies. 

A particularly contentious debate has centred on the in- 
ner slope of the density profile inferred from rotation curves 
of low surface brightness (LSB) disk galaxies. Many au- 
thors have argued for the presence of a constant density 
'core' at the centre of galactic dark matter halos, in direct 
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contradiction with the 'cuspy', centrally divergent density 



profi l e of simulated CDM halos (see, e.g ., Flores fc Primack 
19941 : lMoorelll994l : lMcGaugh fc de B lok 1998; de Blo k et al 



20011 ) . This conclusion relies on the assumption that the 



disk is in circular motion and, therefore, that the observed 
rotation velocity is a fair tracer of the circular velocity of 
the halo. This assumption holds only for thin, cold, gaseous 
disks in spherically symmetric potentials and, while this is 
a convenient assumption for preliminary studies, it has long 
been known that CDM halos are not spherically symmet- 



Warren et al 



19921; 



ric objects llFrenk et al.l 19881; Dubinski fc Carlberd 



Thomas et al 



19981; iJing fc Sutol 
" 20061 ). 



1991 



2002 



Bailin fc Steinm ctz 2005; Bct t et al 

Recently, lHavashi fc Navarro! (|2006l , hereafter HN) in- 
vestigated the kinematics of a gaseous disk in a triaxial halo 
potential using analytic solutions for closed loop orbits. The 
shape of the disk rotation curve, as measured by a long-slit 
spectrum, depends on the orientation of the slit, the inclina- 
tion of the disk, and the radial dependence of the departures 
from spherical symmetry. HN show that significant devia- 
tions from circular motion may result even in very mildly 
triaxial potentials, since the orbital shape is controlled by 
the ratio of escape to circular velocities, a quantity that in- 
creases steadily inwards in CDM halos. 

Few signatures of the halo triaxial structure — which 
may be recognised in 2D velocity fields — are discernible 
in long-slit observations. Linearly-rising "core-like' rotation 
curves may occur if (i) the slit samples velocities near the 
long axis of the (elliptical) closed orbits and if (ii) the devia- 
tions from spherical symmetry in the halo potential become 
more pronounced towards the centre. As noted above, the 
shapes of CDM halos have been studied previously; however, 
these studies have focused on inertia-tensor shapes rather 
than on the potential, thus it is still unclear whether the 
radial dependence of halo triaxiality is quantitatively con- 
sistent with the latter requirement. 

Establishing the typical radial dependence of triaxial- 
ity in the potential of CDM halos is therefore one of the 
main concerns of this paper. The interest of this study goes 
beyond the topic of LSB rotation curves, since the non- 
spherical nature of CDM halos affects the interpretation of 
a variety of observati onal probes, includi ng the kinematics 
of polar ring galaxies ( |Sackett _ et al.l[l 994|), the warping and 
flarin g of galactic disks (jTeuben 199 it Oiling & Mcrrificld 
2000), and the morphology of tidal streams. In a spherical 
potential, tidal streams remain confined to a plane, whereas 
the pl ane of motion will precess gradually in non-spherical 
halos. Ilbata et al.l i|200ll ) use this to interpret tidal debris 
from the Sagittarius dwarf galaxy, and conclude that the 
halo of the Milky Way is spherical and ther efore possibly 
in conflict with CDM predictions. However, iHelmil (12004 ) 
argues that the Sagittarius stream may be too dynamically 
young to be sensitive to the shape of the dark matter halo, 
and hence may be consistent with a Milky Way halo as tri- 
axial as a typical simulate d CDM halo (see also lLaw et al.l 
120051 ; I Johnston et al.ll2005h . 

In galaxy clusters, the hot intracluster gas is in equi- 
librium with the gravitational potential of the cluster, and 
the isodensity surfaces of the gas coincide with the isopoten- 
tial surfaces of the cluster halo (????). Observations of X- 
ray emission and Sunyaev-Zel'dovich (S-Z) decrement mea- 
sure the integrated density along the line-of-sight though the 



cluster, therefore the shape of the X-ray isophotes and S-Z 
contours reflects the shape and orientation of the underlying 
dark matter halo potential. Models of the projected surface 
brightness profiles as a function of cluster orientation and 
triaxiality have been developed by ? and ?, but these are 
based on a constant triaxiality in the gas isodensity sur- 
faces, which may not be an accurate representation of the 
distribution of gas in equilibrium with a realistic CDM halo 
potential 

There is therefore considerable interest in firming up 
the theoretical expectation for the shape and radial depen- 
dence of the gravitational potential in CDM halos, and this 
paper aims to provide guidance for such modelling. We fo- 
cus here on the shape of the isopotential surfaces of galaxy- 
sized dark matter halos using high-resolution cosmological 
N-body simulations. 

The outline of the paper is as follows. In Section [5] we 
describe the simulations and methods we use to calculate 
halo shapes. In Section [3] we present our measurements of 
the shapes of halos and of the internal alignment of halo 
shapes. In Section|4]we investigate the flattening of the mass 
distribution required to reproduce the radial dependence of 
the halo potential shapes, and we present a simple fitting 
formula to approximate the potential of a triaxial CDM halo. 
Finally, we apply this result in Section [S]to the kinematics of 
disks in realistic triaxial CDM halo potentials. We conclude 
with a brief summary in Section [6] 



2 SIMULATIONS AND METHODS 

This study is based on a suite of cosmological N-body simu- 
lations of seven Milky Way-sized galaxy halos and four dwarf 
galaxy-sized halos. The concordance ACDM cosmology is 
adopted for these simulations, with flo = 0.3, £Ia = 0.7, 
o"8 = 0.9. and either h — 0.65 (runs labelled Gl, G2 and 
G3) or h — 0.7 (the rest). These s imulations make u se of 
the techniques descri bed in detail in lPower et al.l (|2003l ) and 
iNavarro et al.l l|2004l ) for resimulating halos at high resolu- 
tion using nested regions of particles with different mass 
resolution and gravitational softening lengths in order to 
maximize the number of particles which end up in the tar- 
get halo at z — 0. 

Each halo contains about N ~ 10 6 particles within the 
virial radius, V2oo, defined as the radius of a sphere of mean 
density equal to 200 times the critical density for closure. 
The spherically-averaged mass profile of these halos is ro- 
bust down to r conv ~ 0.01 V2oo accordin g to the convergence 
criteria described in iPower et al] (|2003l ). Table [TJ summaries 
the main properties of the simulated halos; a full description 
of the numerical details of these simul ations and an analysis 
of t he halo mass profiles i s presented in lHavashi et al.l (|2004| ) 
and lNavarro et all l|2004h . 

Figure [T] shows one galaxy-sized halo from our set of 
simulations. The upper panels show the halo with its parti- 
cles coloured according to the values of their local density 
(left) and gravitational potential (right). The local density 
of each particle is computed by averaging over its 64 nearest 
neighbour particles with a spline kernel similar to that used 
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in Smoothed Particle Hydrodynamics (SPH) calculations^] 
The gravitational potential of each particle is computed us- 
ing all particles within 4 r^oo of the halo centre, determined 
using an iterative technique in which the centre of mass of 
particles within a shrinking sphere is c omputed recursivel y 
until a few thousand particles are left (|Power et al.ll2003h . 
The bottom panels plot the values of the density and gravi- 
tational potential of each particle as a function of its distance 
from the halo centre. 

The most striking feature of this figure is that the 
wealth of halo substructure readily apparent in local density, 
is much less prevalent in the gravitational potential. This is 
not only because the potential is an integral (and therefore 
more uniform) quantity, but also because the total amount of 
mass associated with substructure is typically only 5 — 10% 
of the mass within r2oo (see, e.g., iDe Lucia et al.l l2004h . 
The lower left panel of Figure [T] shows that particles within 
substructures (subhalos) appear as sharply defined "spikes" 
which can exceed the mean local density at the radius of the 
subhalo by many orders of magnitude. This can cause un- 
desirable bia s es wh en measuring halo shapes. For instance, 
Ijing fc Sutol (|2002) estimate shapes using isodensity sur- 
faces, a procedure that requires a careful, and somewhat 
arbitrary, removal of substructure as a preliminary step, and 
that may thus introduce ambiguities in the results. 

In comparison, as emphasised bv lSprineel et alJ l|2004h . 
the isopotential surfaces are much smoother and relatively 
insensitive to the presence of substructure. The lower right 
panel of Figure [1] shows that the gravitational potential of 
particles in subhalos deviates from the mean value at the 
subhalo radius by at most a factor of two. Furthermore, it 
is the gravitational potential and not the local density that 
is the more relevant quantity in most dynamical studies of 
halo structure. 

In order to measure the three-dimensional shape of a 
ha lo's isopotentia l surfa ces, we adopt the method described 
bv lSpringel et all (|2004h . The gravitational potential is first 
computed on three uniform grids covering orthogonal planes 
which intersect at the location of the minimum gravitational 
potential of the halo. In our standard set-up, these meshes 
have 1024 2 cells and extend to a distance of 2 T2oo- We use 
a hierarchical tree algorithm to compute the potential effi- 
ciently and include all the mass out to a radius of 4 V2oo when 
processing an individual halo, while more distant particles 
are ignored. 

If the isopotential surfaces are ellipsoids, then their in- 
tersections with the three principal planes are ellipses, and 
the full shape information of each ellipsoid can be recovered 
uniquely from the shapes of the intersections. Compared to 
a full 3D grid, using just three planes has the important ad- 
vantage of allowing a much finer mesh while reducing the 
memory and computational requirements by several orders 
of magnitudes. 

In our actual fitting procedure, we first determine in 
each plane the intersections of 360 radial rays with constant 
angular separation with the isopotential ellipses for a chosen 
value of the potential. We define the potential at every point 
in the planes by bi-linear interpolation in the corresponding 
two-dimensional grid. We adopt the centre-of-mass of the 

1 See http: //www-hpcc . astro . Washington. edu/tools/smooth. html 



resulting set of points as the centre of the isopotential ellip- 
soid, and identify its principal axes with the eigenvectors of 
the moment-of-inertia tensor of the sample of points. To de- 
termine the axis lengths (a, b, c) , we first express the coordi- 
nates of the points in the orthonormal basis of the principal 
axis relative to the ellipsoid's centre. Denoting these coordi- 
nates with (xi,yi,Zi) we then define a normalized radius r, 
for each point as 

7\ ~~d 2 w 1?' 1 ' 

If a point lies on the ellipsoid, r; would just equal the Carte- 
sian distance \J x'\ + yf + zf, and the left-hand-side of equa- 
tion fl]) would be unity. To determine the best-fitting axis 
lengths, we therefore minimize the quantity 

s = J2( r >-^ x *+y* + z ?) 2 ( 2 ) 

with respect to (a,b,c). In practice, we use a Newton- 
Raphson method to alternatingly find the roots in §^ =0, 
= 0, and ^7 = 0, cycling through the equations until no 
further reduction of S can be achieved. 

Figure [2] shows the results of this procedure applied to 
galaxy halo G7. The right panels show the ratio of the semi- 
major axes of the best fitting ellipsoids, b/a and c/a, where 
c < b < a, as a function of the radius r' — \J a 2 + b 2 + c 2 . 
The radius is plotted in units of the NFW scale radius, r s 
(where the logarithmic slope of the density profile is equal 
to the isothermal value of —2), determined by fitting the 
spherically-averaged density profile with the NFW fitting 
formula. 

The top panels show the highest resolution realization, 
G7/256 3 , which has N200 — 3.5 x 10 6 particles within r200 
(shown by a circle in the left panels). Both the minor-to- 
major (c/a) and intermediate-to-major (b/a) axial ratios de- 
crease gradually inwards, with evidence of a sharp decrease 
inside the scale radius. The top right panel also illustrates 
the uncertainty in the shape measurements by plotting the 
axial ratios calculated using two sets of orthogonal planes 
rotated by an arbitrary angle. The axial ratios measured in 
both cases are in good agreement and differ by <C 3% at 
all radii, confirming the applicability of approximating the 
isopotential surfaces as ellipsoids. 

The middle and lower panels of Figure [2] show the re- 
sults for two lower resolution simulations of the same halo 
G7. The realizations labelled N = 128 3 and N = 64 3 have 
N200 — 3.3 x 10 5 and 4.2 x 10 4 , respectively. The axial ratios 
of G7/128 3 are in good agreement with those of the high res- 
olution N = 256 3 run, with maximum deviations of ~ 4%. 
In the lowest resolution run, however, although the general 
trend in the axial ratios is reproduced, large deviations are 
apparent at small radii. Near the centre, G7/64 3 is signif- 
icantly more spherical than higher resolution realizations, 
suggesting that the steep inward increase in asphericity seen 
in G7/128 3 and G7/256 3 is a result of the inner structure 
of the halo, which is poorly resolved in the case of G7/64 3 . 
The two-body relaxation timescale depends on the number 
of particles and the gravitational softening length and there- 
fore is shortest near the centre. This may result in a more 
spherical distribution of particles, which in turn can cause 
the potential to be more spherical near the centre of low 
resolution halos. We investigate this possibility in Section [4j 
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but for the rest of the analysis we use only simulations of 
resolution comparable to the 256 3 realization of halo G7. 



3 HALO SHAPES AND INTERNAL 
ALIGNMENT 

3.1 Radial dependence of isopotential shapes 

Because the monopole dominates in the outer regions of 
a centrally concentrated mass distribution, the isopotential 
surfaces are expected to become more spherical in the outer 
regions. The precise radial variation of the shape depends on 
the distribution of mass within the halo and on the radial 
dependence of its flattening. 

Figure [3] shows the axial ratio profiles for all of the sim- 
ulated halos in our sample along with the average profile 
(computed after scaling each profile to its scale radius, r 3 ) 
with la error bars to indicate the scatter. The slow approach 
to spherical symmetry in the outer regions is clear in all 
cases (the average minor-to major axial ratio at r = 5 r s is 
c/a ~ 0.9), as is the sharp inward increase in asphericity in- 
side the scale radius. Near the centre, c/a and b/a approach 
similar values, 0.72 and 0.78, respectively, but the system 
becomes more spherical and less axisymmetric in the outer 
regions. Figure [3] also shows the triaxiality parameter, de- 
fined as T = (a 2 - b 2 )/(a 2 - c 2 ), where T = (1) for a 
perfectly oblate (prolate) spheroid. On average, T > 0.5 for 
r < r s indicating that halos tend toward prolate shapes near 
the centre. The c/b profiles remain relatively flat from the 
centre outwards, increasing from a central value of ~ 0.92 
to ~ 0.94 in the outer parts of the halo. 

A more detailed view of this inner region is shown in 
Figure|3J where we plot b/a versus c/b at four different radii, 
r = 0.1,0.25,0.5, and 1.0 r s . Note that a perfectly oblate 
spheroid corresponds to b/a = 1, whereas a prolate spheroid 
corresponds to c/b = 1. There is a clear trend for halos 
to become increasingly prolate near the centre. At r — r s , 
five out of eleven halos are more prolate than oblate, but at 
r = 0.1 r s , all but two of the halos are more prolate than 
oblate. 

The significant variation in the shape of the halo po- 
tential with radius may complicate detailed comparisons of 
the shapes of simulated halos with those inferred from X-ray 
and S-Z observations. In the model of ?, for example, this 
variation is neglected in order to simplify the calculation of 
the projected surface brightness profile of a galaxy cluster. 
The large range in axial ratio values at r <C 2 r s suggests 
that comparing axial ratios measured at large radii might 
be a more straightforward test of the distribution of shapes. 
Observing an increase in the ellipticity of X-ray isophotes or 
S-Z isodecrement contours towards the centre in individual 
cluster systems would also provide strong evidence for the 
existence of an underlying CDM halo potential. 



3.2 Alignments 

The alignment of the isopotential surfaces as a function of 
radius is illustrated in Figure [5] which shows the cosine of 
the angle between the principal axes of the isopotential at 
the "converged" radius and the corresponding axis at larger 
radii. Since orientations between some axes will be poorly 



determined in cases where b ~ a or c ~ b, we plot the results 
using solid (dotted) curves in Figure [5] when the length of 
the axis plotted differs by more (less) than 5% from the 
lengths of the other axes. Taking this into account, we find 
that in most cases where the axes are well determined the 
isopotential surfaces appear to be reasonably well aligned 
as a function of radius, with the possible exception of halos 
G2, G3, and D3. 

The bottom- right panel of Figure [S] also shows the rel- 
ative alignment between the potential and the angular mo- 
mentum of the halo, by plotting the cosine of the angle be- 
tween the minor axis and the total angular momentum, com- 
puted using all particles within T2oo- We find that the minor 
axis tends to be aligned with the angular momentum vector 
in most halos; in six of the eleven halos, the alignment be- 
tween the minor axis and the angular momentum vector is 
better than 25° out to ~ 5 r„. 

The latter result has important implications for the dy- 
namics of disk galaxies since baryons and dark matter are 
expected to have acquired similar (specific) angular mo- 
menta, and, therefore, the rotation axis of a galactic disk 
would be aligned with the halo's net angular momentum 
fsee. e.g..lFall fc Efstathioulll98rj|;lyan den Bosch et al.ll2002l ; 



iNavarro et alj|2004 lAbadi et alj|2006h . Such a disk would 
be confined to a plane containing the major and intermedi- 
ate axes of the halo. In a nearly prolate halo the potential 
in such plane would deviate significantly from circular sym- 
metry. 

We show this explicitly in Figure [6] where we plot the 
potential of halo D3 on the plane containing its major and 
intermediate axes. The isopotential contours are clearly non- 
circular but are well-approximated by ellipses. The four pan- 
els in Figure [6] show the halo on increasingly small scales in 
order to highlight the increase in the ellipticity of the isopo- 
tential contours towards the centre of the halo. We investi- 
gate the consequences of this for the dynamics of a gaseous 
disk in Section [S] 



4 A SIMPLE MODEL FOR HALO 
TRIAXIALITY 

The radial dependence of the shape of the potential is dic- 
tated by both the flattening and the mass profile of the halo. 
Here we explore the effects of each in order to shed light on 
the steep radial dependence of the halo asphericity discussed 
in the previous section. 

To this aim, we generate a spherically symmetric NFW 
halo model with N = 10 6 particles out to a maximum ra- 
dius of 10 r s . We set the gravitational softening length to 
0.05 r s . We turn this into a prolate halo of constant flatten- 
ing by multiplying the x- and y-coordinates of all particles 
by a factor of 0.4. We then use the ellipsoid fitting method 
described in the previous section to measure the shapes of 
the isopotential surfaces as a function of radius. We also 
measure the shape of the mass distribution by diagonaliz- 
in g the i nertia tensor according to the iterative procedure 
of lKatzl 119911 ). This method measures the shape using all 
particles within an ellipsoid whose volume is approximately 
equal to the volume of a sphere of radius r. 

Figure [7] compares the shape of the self-similar pro- 
late NFW halo (left panel) and halo G6 (right panel). The 
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dashed curve (without symbols) in the left panel shows the 
axial ratio measured using the inertia tensor; as expected, 
a constant value of ~ 0.4 is recovered for the prolate NFW 
model. The potential axial ratio (dotted curve without sym- 
bols) is not constant, but the radial dependence is much 
gentler than seen in halo G6. Within r s the shape of the 
potential of the prolate NFW model remains more or less 
constant at ~ 0.65; and becomes gradually more spherical 
outside r s . This differs from the radial behaviour of the ax- 
ial ratios in halo G6 (right panel) where both the poten- 
tial and inertia shapes within r 3 increase from the centre 
outwards. These results indicate that the sharp increase in 
the asphericity of the halo potential is inconsistent with a 
constant flattening in the mass distribution and reflects the 
presence of a radially varying distortion in the mass distri- 
bution as well. 

A convenient approximation of the shape profiles is pro- 
vided by the following formula: 



'b c 
log — or — 

a a 



tanh 7 log — 



1 



(3) 



Here a parameterizes the central value of the ratio, such that 
(b/a)o or (c/a)o is given by 10 -2 a , r a indicates the char- 
acteristic radius at which the axial ratio rises a significant 
amount from its central value, and 7 regulates the sharpness 
of the transition. This fitting formula gives a reasonably ac- 
curate description of the axial ratio profiles of both the mass 
distribution and the isopotential surfaces in the simulations, 
as shown by the fits to the G6 shape profiles in Figurc[7](lines 
through symbols in its right panel). The best fit parameter 
values for fits to the mass and isopotential shapes are given 
in Tabled] 

An NFW halo model with a radially varying flattening 
is a much better approximation to the radial behaviour of 
the shapes seen in the simulations. This is shown in the left 
panel of Figure [7] where the same (initially spherical) NFW 
model is modified by multiplying the x- and y- coordinates 
of particles at radius r by eq. (J3J) evaluated with (a, r a ,"f) 
= (0.2, 1.0 r s , 1.3). This choice is motivated by the best fit 
to the b/a profile of the mass distribution of halo G6. 

The axial ratio profiles of the resulting halo model are 
shown by the symbols in the left panel of Figure [7] The dis- 
tortion in the shape of the isopotential surfaces of this halo 
is clearly reminiscent of that of halo G6. The agreement is 
not exact, however, because the inertia tensor measures the 
shape using all particles within a given volume and there- 
fore does not correspond to a transformation applied to the 
coordinates of particles as a function of radius. However, the 
model does demonstrate that a distortion in the mass dis- 
tribution which increases towards the centre of the halo can 
reproduce the radial variation observed in the isopotential 
shapes of halos like G6. 

We also find that in all cases, the potential axial ratio 
profile becomes flat, or increases slightly, toward the centre 
of the halo. In G6 and the radially varying NFW model, this 
coincides with a similar trend in the mass axial ratio profile. 
However, in the NFW model with a constant flatenning, the 
shape of the potential becomes more spherical near the cen- 
tre, but the shape of the mass distribution remains constant 
at all radii. This suggests that the finite gravitational soften- 
ing length may be responsible for making the potential more 
spherical near the centre, although two-body relaxataion ef- 



fects may also play a role near the centres of simulated halos 
like G6 and the low resolution versions of G7 shown in Fig- 
ured 

We note that the ratio of the potential axial ratio to the 
mass axial ratio varies as a function of radius for both the 
constant and radially varying flatenning models. In terms 
of the eccentricity parameter used in the models of ? and 
?, e — wl — (b/a) 2 ), the ratio of potential eccentricity to 
the mass eccentricity varies from ~ 0.8 to ~ 0.6, and ~ 0.8 
to ~ 0.3, in the constant and radially varying cases, respec- 
tively. In comparison, the eccentricity ratio of the ? model 
varies from 0.77 to 0.4, however, ? approximate this with 
a constant value of 0.49 in order to calculate the projected 
surface brightness distribution of their model. Since this ra- 
tio varies over roughly the same range for the model with 
a radially varying flatenning, we conclude that the assump- 
tion of a constant ratio is no worse in this case than it is for 
models with constant flatenning. 

The best fit parameters using eq. <(3j to approximate the 
potential axial ratio profiles of all our simulated halos are 
listed in in Table [TJ and shown in Figure [5] Although there is 
substantial halo-to-halo scatter, several trends are apparent. 
The large symbols with "error bars" in Figure [8] show the 
average value of each parameter and the 1-sigma scatter. 
The central axial ratios scatter about (6/a)o = 0.78 ± 0.08 
and (c/o)o = 0.72 ± 0.04. The average transition radius, 
expressed in units of the scale radius, r a /r s , is 1.2 ± 1.2 
for the b/a profiles (excluding halo G5 which has an almost 
flat b/a profile, and therefore an extremely large and poorly 
defined value of r a ) and 3.0 ± 2.4 for the c/a profiles. The 
average value of 7 is 1.4±0.8 for the b/a profiles and 1.1±0.3 
for the c/a profiles. On average, c/a profiles tend to increase 
more gradually than b/a profiles, and profiles become more 
spherical and less axisymmetric at larger radii. 

Our results may be used to construct a model for the 
gravitational potential of a triaxial CDM halo by using 
eq. (|3} to modify the potential of a spherical NFW model. 



$nfwW 



r{r) 



— + 4- 

a(r) J \b(r 
GM200 ln(l + r/r 
r s /(c2oo) r/r s 



+ 



z 
c(r) 



(4) 
(5) 
(6) 



where the concentration parameter C200 = rwo/r s and 
f(u) — ln(l + u) — u/(l + u). Since eq. Q describes only the 
ratios between the lengths of the principal axes, the normal- 
isation is set by assuming that at each radius r, the volume 
of the ellipsoid with semi-axis lengths a(r),b(r) and c(r) is 
equal to 47rr 3 /3, i.e., a(r)b(r)c(r) = r 3 , which implies 



a(r) =r (b/a)- 1/3 (c/a)- 1/3 , 



(7) 



where b/a and c/a are functions of r given by eq [3] This 
ensures that the spherically-averaged potential of the triax- 
ial model will be similar to that of a spherically symmetric 
NFW model with the same mass and scale radius. 

We show that this model provides a good description 
of the potential of simulated halos in Figure [9] where we 
compare the gravitational potential of particles in halo G5 
plotted against their radius, r, and against the reduced ra- 
dius, r , given by eq. ([5]). The reduced radius is calculated 
by rotating the halo so that its principal axes are aligned 
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with the x-, y-, and z-coordinate axes. The bottom panels 
of this figure show that the dispersion in the potential pro- 
file is significantly reduced when expressed in terms of the 
reduced radius; in the inner regions, r < 2 r„, the dispersion 
is reduced by a factor of two. At larger radii, the presence 
of substructure dominates the scatter in the potential and 
the triaxial modelling has less impact. 

We conclude that taking into account halo triaxiality 
provides a substantial improvement over spherically sym- 
metric CDM halo models. The fitting formula given by 
eq. ((3} along with parameter values in the ranges we find 
for simulated halos, listed in Table [T] provides a simple ana- 
lytic model for a realistic, triaxial CDM halo potential. This 
model may be of prac tical use for N-body simu lations of tida l 
streams like those of llbata etafl (|200lh and iHelmil i|2004h . 
and in calculations of the X-ray and S-Z surface brightness 
distributions of triaxial galaxy clusters as in ?. 



5 DISK KINEMATICS IN TRIAXIAL CDM 
HALOS 

We explore now the kinematics of disks in CDM halos with 
radially varying triaxiality. We assume that the disk lies in 
one of the principal planes of the potential, and therefore 
the problem is reduced to two dimensions. Our approach 
follows closely the formalism presented in HN, where ana- 
lytic closed loop orbit solutions are found to describe the 
orbits in a thin, filled, gaseous disk. These solutions use the 
epicyclic approximation to find closed orbits in a potential 
of the form: 



decreasing nearer the origin. This behaviour can be approx- 
imated by the following fitting formula: 



$o(i?) + & m (R) cos(m0 o ) 



(8) 



where &o(R) is the unperturbed (circularly-symmetric) po- 
tential, <!> m (R) is a stationary perturbation to that potential, 
and 4> is the azimuthal angle. For m = 2 and small pertur- 
bations, the perturbation is periodic in 180° and the isopo- 
tential contours are roughly elliptical in shape. The closed 
loop orbits solutions are also roughly elliptical, but are ori- 
ented with their major axes perpendicular to those of the 
isopotential contours. 

The tangential velocity along the orbit oscillates about 
the circular velocity of the unperturbed potential, V C (R) = 
(Rd&o/dR) 1 / 2 , and the deviations from V c will be max- 
imal at pericentre and apocentre of the orbit. Long-slits 
aligned with these directions would yield rotation curves 
whose shape will systematically deviate from the true circu- 
lar velocity profile. 

HN showed that even small perturbations to the po- 
tential (i.e., $ m <C <E>o) may result in substantial deviations 
to the orbital shapes and thus to the tangential velocities; 
rotation curves that differ from the halo's V c profile may 
thus be matched by choosing a specially tailored function 
& m (R)- In particular, HN compute the perturbation needed 
to match a pseudo-isothermal velocity profile (correspond- 
ing to a isothermal sphere with a constant density core) in 
a cuspy NFW profile. 

This is shown by the dashed curve in the upper right 
panel of Figure 1101 The perturbation has a characteristic 
shape: it increases inwards to a well defined peak that oc- 
curs well inside the scale radius of the NFW halo, before 



<F, 



R 



exp 



R 



0) 



where r„ 



0.098 r s and a„ 



9.8 x lO -3 . Note that 



r m <C r a , as required to affect the rising part of the NFW 
rotation curve (which peaks at ~ 2r s ), and that the over- 
all perturbation is relatively minor (everywhere less than 
0.4%). 

Is this perturbation consistent with the radially- varying 
triaxial structure of CDM halos described in the previous 
section? We calculate <I> m for a simulated halo with the fol- 
lowing procedure. We compute the potential on concentric 
rings on the plane containing the intermediate and major 
axes (the most likely one to contain the disk given the align- 
ment between angular momentum and minor axis) . The top 
left panel of Figure [10] shows the potential plotted versus 
azimuthal angle for halo G4 on a ring of radius 0.5 r s . The 
potential is periodic in ir radius, and is reasonably well fit 
by a sinusoid with three free parameters which describe its 
phase, amplitude and mean value, i.e., Acos(2(/> + <f>o) + B. 
At each radius, we fit such a sinusoid to the potential and 
use the best fit parameters to calculate the magnitude of the 
perturbation, |$ m /^>NFw| — \ A/B\. The result is shown for 
halo G4 as the dot-dashed in the top right panel of Figure fTOl 

The magnitude of the perturbation peaks at r ~ 0.2 r 3 
and is about 2.5 times as large as the minimum needed to 
reconcile a "core-like" rotation curve with a cuspy NFW pro- 
file. Unlike the solution presented in HN, however, the mag- 
nitude of the perturbation for halo G4 does not decrease to 
zero at large radii and in this case the closed orbit solutions 
based on the epicyclic approximation break down. We ap- 
proximate the perturbation of halo G4 by fitting it over the 
range R < 0.5 r s with eq. (|9}. This fit is shown as the solid 
line labelled $fl t in the top right panel of Figure 1101 and the 
fit parameter value 0.027 and r m = 0.21 r s . 

The bottom left panel of Figure \W\ shows the b/a ax- 
ial ratio profiles corresponding to halo G4 and the per- 
turbed NFW potential models. The axial ratio profile of 
the f&fit perturbation is quite similar to that of halo G4 
and is well fit by eq. © with parameter values (a, r a ,j) = 
(0.079,0.24 7-3,2.68), compared to (0.062,0.22 r s , 1.85) for 
halo G4. 

The lower right panel of Figure [TU] shows as open cir- 
cles the rotation curve that would be obtained by a long 
slit sampling tangential velocities along the major axis of 
the disk (cj> a — 0° in the notation of HN) in the per- 
turbed potential given by the fit to halo G4. The slope 
of the inner rotation curve is significantly modified from 
that of the initial NFW profile, and is well fit by a pseudo- 
isothermal velocity profile (solid line), given by V^ (r) = 
V^[l — (r COIC /r) tan _1 (r-/r corc )], where Vac is the asymptotic 
velocity and r corc is the radius of the constant density core 
in this model. The best fitting value of the core radius in 
this case is r corc = 0.38 r s , comparable to the characteristic 
radius of the perturbation, r m m 0.45 r a « 0.24 r s . 

Extending this comparison to the other halos in our 
simulation suite is not straightforward, because the rather 
large deviations from spherical symmetry present in many 
cases preclude the use of the epicyclic approximation to com- 
pute closed orbits. Our set of 11 simulations is also too small 
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to allow for a statistically meaningful comparison between 
simulated halos and LSB rotation curves. The case of G4 is 
nevertheless encouraging, as it signals that perturbations of 
the right magnitude and the right shape can result in "core- 
like" rotation curves such as those observed in some LSB 
galaxies. 



6 CONCLUSIONS 

We have examined the three-dimensional shape of the grav- 
itational potential in seven Milky Way-sized halos and 
four dwarf galaxy-sized halos simulated in the concordance 
ACDM cosmo l ogy. W e use the ellipsoid fitting method of 
ISpringel et alj (120041 ) to measure the shape of isopotential 
surfaces as a function of distance from the centre of the 
halo. The isopotential surfaces are well described by con- 
centric ellipsoids whose principal axes remain aligned across 
the whole halo, and whose minor axes tend to be parallel to 
the angular momentum of the halo. 

We find that the axial ratios of the isopotential surfaces 
decrease sharply inwards inside the characteristic scale ra- 
dius of the halo mass profile, approaching central values of 
(c/a)o = 0.72 ±0.04 and (b/a) = 0.78 ±0.08. The potential 
becomes more spherical and less axisymmetric from the cen- 
tre outwards, as c/a increases more gradually than b/a. Such 
radial dependence is well captured by a simple fitting for- 
mula, eq. (0, where the characteristic radial scale for c/a is 
about twice that of b/a. Incorporating the radial variation 
in shape in models of CDM halos represents a significant 
improvement over the analytic halo potentials generally im- 
plemented in N-body simulations. 

We use the analytic closed lo op orbit solutions pre- 
sented in lHavashi fc Navarrd (|2006T ) to investigate the kine- 
matics of thin, filled gaseous disks embedded in such halos. 
We find that the radially-dependent ellipticity of the po- 
tential on the principal planes of the halo generally resem- 
bles the perturbation needed to reconcile "core-like" rota- 
tion curves with cuspy halo profiles. The amplitude of the 
ellipticity in the potential is typically larger than needed to 
obtain linearly-rising long-slit rotation curves according to 
the treatment presented in HN. 

This is important, as the simulations do not include 
baryons, and the assembly of a baryonic disk will tend to 
circularize the potential, by (i) inducing ch anges in the ac- 
tual mass distribution of the dark halo (|Dubinskil [l99i 
iGustafsson et alj|2006l . Abadi et al, in preparation), and by 
(ii) "opposing" the halo ellipticity — the disk would itself be 
elliptical but rotated by 90° relative to the halo potential. 
The magnitude of the corrections implied by these two ef- 
fects is at present quite uncertain, and will require simula- 
tions more precisely tailored to reproducing the formation 
of a system with the observed properties of an LSB galaxy. 

Because of these uncertainties it would be premature to 
conclude that halo triaxiality is the main cause of both the 
wide variety in LSB rotation curve shapes and the existence 
of rotation c urves suggestive of the presence of constant- 
density cores (jHavashi et al.ll2004[ ) . Still, our results demon- 
strate that radially-varying ellipsoidal potentials are ex- 
pected to be the rule rather than the exception in CDM 
halos; and indicate that this "natural" feature of CDM halo 
structure needs to be included in more sophisticated analy- 



ses of LSB rotation curves and of other observational probes 
of the shape of CDM halos. 
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Table 1. Main parameters of simulated halos 
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Figure 1. Upper panels: Halo G3/256 3 coloured by local density (left) and potential (right). Color scheme is logarithmic for density 
and linear for the potential. The virial radius, r200i is shown by the dashed circle. Lower panels: Local density (left) and potential (right) 
of particles as a function of radius. Thin solid lines show the average local density (left panel) and best fit NFW potential (right panel). 
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Figure 2. Halo G7 simulated at three different resolutions. The number of particles varies by a factor of 8 between each realization. 
Left panels: Halo particles coloured by gravitational potential. Right panels: Axial ratios b/a and c/a measured by fitting ellipsoids to 
isopotcntial contours on three orthogonal planes through the halo as a function of the elliptical radius r' = (a 2 + b 2 +C 2 ) 1 / 2 . The vertical 
dashed line shows the minimum converged radius in the mass profile, r conv , in each simulated halo. The top right panel shows the axial 
ratios calculated using two sets of orthogonal planes rotated by an arbitrary angle. The axial ratios measured in both cases are in good 
agreement and differ by < 3% at all radii. 
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Figure 3. Isopotential shapes of all halos as a function of radius, where c, 6, and a, are the lengths of the minor, intermediate and 
major axes, respectively, and the radius is given in units of the NFW scale radius of the density profile, r s . The triaxiality parameter 
T = (a 2 — b 2 )/(a 2 — c 2 ) is equal to (1) for a perfectly oblate (prolate) spheroid. Points with error bars show the mean ±1 a of all axial 
ratio profiles. Halos clearly tend to become more aspherical towards the centre, a gradual trend that becomes increasingly steep inside 
r s . On average, halos are very nearly prolate at the centre, but become more spherical and less axisymmetric in the outer regions. 
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Figure 4. Intermediate-to-major (b/a) versus minor-to-intermediate (c/b) axial ratios at four different radii within r s . Perfectly oblate 
(prolate) halos have b/a = 1 (c/b = 1) and "maximally triaxial" halos have b/a = c/b. Halos tend to become more prolate inwards. 
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Figure 5. Alignment of principal axes as a function of radius. Vector a(r') represents the unit vector along the major axis of the halo, 
and do = a(r' ~ r conv ). Dotted lines indicate radii where c/b > 0.95 and/or b/a > 0.95. For these nearly axisymmetric systems two of 
the axes' directions are poorly determined. In eight of the eleven halos, the principal axes are well aligned with radius throughout the 
main body of the halo. In most halos j r 200 ' c(r) ~ 0.9, indicating alignment of 25° or better between the minor axis and the angular 
momentum vector of the halo. 
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Figure 6. Gravitational potential of halo D3 on the plane perpendicular to its minor axis, corresponding to the plane in which a disk 
might form. The radii of the outer dotted circles are 4, 2, 1, and 0.5 r s in the top left, top right, bottom right panels, respectively, where 
r s = 2.94 h~ 1 hpc. The radius of the inner circles is 0.25x that of the outer circle. The ellipticity of the isopotential contours (solid 
curves) clearly increases towards the centre of the halo. 
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Figure 7. Le/t panel: Axial ratios of NFW halo models with constant (lines without symbols) and radially varying (symbols) flattening 
in the mass distribution as a function of radius. Right panel: Axial ratios of halo G6 as a function or radius. Solid curves show fits to the 
axis ratio profiles with eq. J3j. Halo models with constant flattening cannot reproduce the axial ratio profiles of simulated halos like G6, 
which require increasing asphericity towards the centre. 
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Figure 8. Values of central axial ratios, (b/a)o and (c/a)o, and parameters r a and 7 obtained by fitting eq. (|3jl to the b/a and c/a 
axial ratio profiles plotted in Figure [3] Large symbols with error bars show the average and standard deviation for all halos. The average 
transition scale for the c/a profiles is about twice that of the b/a profiles, indicating that the former tend to increase more gradually 
than the latter profiles 
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Figure 9. Left panels: Potential versus radius of particles in halo G5. The solid line shows the profile of a best fit NFW potential. The 
lower panel shows the residuals from the fit, a^, as a function of radius. Right panel: Same as left but plotted versus the reduced radius 
r of the best fitting ellipsoid (see eq. [5j. The dashed line in the lower panel indicates the same residuals, cr$, as in the left panel. The 
residuals are significantly reduced at small radii where substructure does not dominate the scatter. 
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Figure 10. Upper left panel: The potential plotted versus azimuthal angle for halo G4 on a ring of radius 0.5 r s . The potential is well fit 
by a sinusoidal function whose mean value and amplitude reflect the magnitude of the perturbation relative to a spherically-symmetric 
potential. Upper right panel: The magnitude of the "perturbing" potential required to yield "core-like" long-slit rotation curves in NFW 
halos, as presented by Hayashi & Navarro (2006) (dashed curve). The dot-dashed curve corresponds to the perturbation calculated for 
halo G4, as measured in the plane that contains the major and intermediate axes. The solid curve represents a fit to this perturbation 
using eq. J9]l- Lower left panel: Axial ratios of isopotential contours as a function of radius. The dashed, dot-dashed, and solid curves 
correspond to the HN solution, halo G4, and the fit to halo G4, respectively. Lower right panel: Rotation curve of a disk in the perturbed 
potential given by the fit to halo G4, produced when a slit samples velocities near the major axis of the disk (open circles). The dashed 
line shows the circular velocity profile corresponding to the unperturbed, spherically symmetric NFW halo. The solid line shows the best 
fitting pseudo-isothermal profile with a constant density core. 



